APS/123-QED 

Lyapunov Analysis of Homogeneous Isotropic Turbulence 

Nicola de Divitiis 

Department of Mechanics and Aeronautics 
University "La Sapienza", Rome, Italy 
(Dated: August 12, 2009) 

Abstract 

The present work studies the isotropic and homogeneous turbulence for incompressible fluids 
through a specific Lyapunov analysis, assuming that the turbulence is due to the bifurcations 
associated to the velocity field. 

The analysis consists in the study of the mechanism of the energy cascade from large to small 
scales through the Lyapunov analysis of the relative motion between two particles and in the 
calculation of the velocity fluctuation through the Lyapunov analysis of the local deformation and 
the Navier-Stokes equations. 

The analysis provides an explanation for the mechanism of the energy cascade, leads to the closure 
of the von Karman-Howarth equation, and describes the statistics of the velocity difference. 

Several tests and numerical results are presented. 
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I. INTRODUCTION 



This work presents a study of isotropic and homogeneous turbulence for an incompressible 
fluid in an infinite domain. The analysis is mainly motivated by the fact that in turbulence 
the kinematics of the fluid deformation is subjected to bifurcations [l| and exhibits a chaotic 
behavior and huge mixing 2| , resulting to be much more rapid than the fluid state variables. 
This characteristics implies that the accepted kinematical hypothesis for deriving the Navier- 
Stokes equations could require the consideration of very small length scales and times for 
describing the fluid motion 3| and therefore a very large number of degrees of freedom. 
Other peculiar characteristics of the turbulence are the mechanism of the kinetic energy 
cascade, directly related to the relative motion of a pair of fluid particles 

mm 

and 

responsible for the shape of the developed energy spectrum, and the non-gaussian statistics 
of the velocity difference. 

The present analysis assumes that the fluctuations of all the fluid state variables are the 
result of the bifurcations of the velocity field. The evolution in the time of these fluctuations 
is calculated with the Lyapunov analysis of the particle equations of motion. 

The first part of the work deals with the representation of velocity difference between two 
fixed points of the space. This is analyzed with the Lyapunov theory studying the motion of 
the particles crossing the two points. This analysis gives an explanation of the mechanism 
of kinetic energy transfer between length scales and leads to the closure of the von Karman- 
Howarth equation Gj. The obtained expression of the function K{r), which represents the 
inertia forces, is in terms of the longitudinal correlation function and its spatial derivative, 
and satisfies the conservation law which states that the inertia forces only transfer the kinetic 
energy p, |7| ■ 

In the second part, the statistics of the velocity difference is studied through the kine- 
matics of the local deformation and the momentum equations. These momentum equations 
are expressed with respect to the referential coordinates which coincide with the material 
coordinates for a given fluid conflguration [s^, whereas the kinematics of the local deforma- 
tion is analyzed with the Lyapunov theory. The choice of the referential coordinates allows 
the velocity fluctuations to be analytically expressed in terms of the Lyapunov exponent of 
the local fluid deformation. The statistics of velocity difference is studied with the Fourier 
analysis of the velocity fluctuations, and an analytical expression for the velocity difference 
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and for its PDF is obtained in case of isotropic turbulence. This expression incorporates an 
unknown function, related to the skewness, which is identified through the obtained expres- 
sion of K{r). The velocity difference also requires the knowledge of the critical Reynolds 
number whose estimation is made in the Appendix B, where the order of magnitude is 
roughly determined through a qualitative analysis of the bifurcations of the velocity field. 

Finally, the several results obtained with this analysis are compared with the data existing 
in the literature, indicating that the proposed analysis can adequately describe the various 
properties of the fully developed turbulence. 

The present analysis only studies the possibility to obtain the fully developed 
homogeneous-isotropic turbulence in a given condition and does not analyze the intermediate 
stages of the turbulence. 

II. LYAPUNOV ANALYSIS OF THE RELATIVE MOTION: CLOSURE OF THE 
VON KARMAN-HOWARTH EQUATION 

In order to investigate the mechanism of the energy cascade, the properties of the relative 
equations of motion between fluid particles are here studied with the Lyapunov analysis. To 
this purpose, consider two fixed points of the space, X and X' (see Fig. [1]) whose distance 
is r = |X' — X| and the motion of two fluid particles which at a given time to, cross through 
X and X'. The equations of motion of these particles are 

^ = u(x,t), ^ = u(x',t) (1) 

At to, a toroidal volume S(to) is chosen which contains X and X', whose geometry and 
position change according to the fluid motion. In Fig. [H 5*^ ^ and R are, respectively, 
the poloidal surface and the toroidal dimension of S which vary with time to preserve the 
volume. The velocity difference associated to X and X' is Au = u(X', t) — u(X, t) and its 
components Am„ = — m„ and Aur = u'^ — Ur, lay on Sp and are normal and parallel 
to r, respectively, whereas Ub is the average of the velocity components along the direction 
normal to Sp. According to the theory [8], for t > to, the trajectories of the two particles 
are enclosed into For sake of simplicity and without loss of generality, we assume that 

R increases with time jsj. The Lyapunov analysis of Eqs. ([T]) leads to 

R ^ R{to) e^(*"*°) (2) 



FIG. 1: Scheme of the relative motion of two fluid particles 



These variations of R are caused of the bifurcations [8| of Eqs. ([T]). Since R rises with time, 
A(r) > identifies the maximal finite scale Lyapunov exponent associated to Eqs. ([1]). 

The equations of motion for ^(t) preserve the volume and can be expressed in terms of 
the velocity components calculated at X and X' j^. These are 

^ {S,R) = (3) 



(4) 



In line with Lamb Q], Eqs and (jl]) represent, respectively, the continuity equation and 
the momentum equations which can be derived from the integral equations of balance over 
S. Into Eq.Q, u is the kinematic viscosity, /3 = O(l)>0isa proper constant, and J is 
related to the time derivative of the kinetic energy and to the viscosity [9|. J is equal to 
zero when u = 0. 
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Equations ([3]) and (jl]) are written in terms of the fluid properties calculated at X and 
X', thus are referred to the Eulerian description of motion 3|, l9|. Substituting Eq. ([2]) into 
Eqs. ([3]) and (jl]), one obtains 



dt 



dt 



\{Au^y + J (5) 



Since A > 0, — 0. Equations describe fluctuations of velocity difference caused by 
bifurcation of Eqs. ([1]) and hold as long as X and X' are both enclosed into This 
condition is verifled if t — to does not exceed very much the order of magnitude of 1/ A ^ . 

In order to obtain the closure of the von Karman-Howarth equation, Eq. ([5]) enter the 
computation of the average of the physical quantity T: 

d , d , 1 dT,, , 1 dTl^ 

-(TrO--MJ--^..--^.. (6) 

where Tij = —pSij + up {dui/dxj + duj/dxi) is the stress tensor. The repeated indexes denote 
the summation with respect to the same indexes, which are i = r,n,b and j = r, n, b. 

According to von Karman gI, T expresses that part of the inertia forces, responsible 
for the transferring of the kinetic energy between the several fluid regions, whose average 
only depends on the current value of the average kinetic energy. In the von Karman- 
Howarth equation, the function K{r) is the average of T. The average is calculated on 
all the pairs of particles which cross through X and X' at the same time. Specifically, 
K{r) is determined substituting Eqs. ([5]) into Eq. IQ, assuming the homogeneity and the 
isotropy and taking into account that (T) does not depend neither on d{uiUi) /dt, nor on 
{dTij/dxju'i + dTlj / dx'jUi) jo, 7|. This immediately identifies 

K{r) = (T) = A u' {g - f) 

(7) 
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where = (uiUi) /3 and / and g are longitudinal and lateral velocity correlation functions. 
Due to the fluid incompressibility, / and g are related each other through g = f + l/2df/dr 
(see Eq. (IS]), Appendix A), leading to the expression 

K{r) = ^u'^ A(r)r (8) 

Considering that K{r) does not directly depend on the viscosity, this expression can be 
also obtained at u = 0. In this case (J) = 0, d{uiUi)/dt = 0, {dTij/dxjU^ + dT-j / dx'jUi) = 
0, and Eqs. ([7]) and ([8]) are again recovered. 

Equation ([8]) states that, the fluid incompressibility, expressed hj g — f 0, represents a 
sufficient condition to state that K{r) ^ 0. This latter is determined as soon as A is known. 
To calculate A, it is convenient to express Au = u(X',t) — u(X, t) in the Lyapunov basis of 



orthonormal vectors E = {ei,£2,e3) associated to Eqs. ([T]) 10, lUl]- The velocity difference 



expressed in E, Av = {v[ — Vi, v'2 — V2, ^3 — fa), satisfies the followiiig equations, which 



hold for times whose order of magnitude do not exceed very much 1/A 8l. Il2| 



f • -Vi = Xi fi, i = 1, 2, 3 (9) 

where fj, Vi and v'^ are, respectively, the components of X' — X, u(X, t) and u(X', t) written 
in E. Then, r and Aur can be expressed in terms of r and Av as 

r = ^ ■ Qf, Aur = ^- QAv (10) 

Into Eqs. ( |T0|) . Q = {{eij)) is the rotation matrix transformation from E to 3?, where Eij is 
the component of £j along the coordinate direction i on 3?, and ^ = (X' — X)/|X' — X|. 

The standard deviation of Am,, is calculated from Eqs. (ITO!) . taking into account that 
Av ^ Af and that Q is fluctuating depending on the pair paths: 

{(Aurr) = (11) 

i=l j = l p=l q=l 

Since A is calculated as the average of the velocity increment per unit distance, it is constant 



with respect the statistics of e^j and Epg 



ll|, thus {X^SijEpq) = X^{eijepq). Furthermore, due 



to isotropy, the Lyapunov vectors fluctuate in such a way that {SijEpq) = As the 

result, the standard deviation of the longitudinal velocity difference is 

{AurY) = (12) 



This standard deviation can be also expressed through the longitudinal correlation function 
/ 

((An,,)') = 2u\l - /(r)) (13) 

being u the standard deviation of the longitudinal velocity. The maximal Lyapunov exponent 
is calculated in function of /, from Eqs. f[T^ and f[T^ 

A(r) = ^V2(l-/(r)) (14) 

Hence, substituting Eq. (031) into Eq. (IHl), one obtains the expression of K{r) in terms of 
the longitudinal correlation function 

Kir) = f (15) 

Thanks to the isotropy, K{r) is a function of r alone. 

Equation f|T5|) represents the proposed closure of the von Karman-Howarth equation. 
This is the consequence of the fact that, the kinetic energy, initially enclosed into S(to), at 
the end of the fluctuation is contained into S(t) whose dimensions are changed with respect 
to S(to) in agreement with the Lyapunov theory. This corresponds to a mechanism of the 
kinetic energy transferring between diverse regions of space which preserves the average 
values of the momentum and of the kinetic energy. Specifically, the analytical structure of 
Eq. ffTSl) states that this mechanism consists of a flow of the kinetic energy from large to 
small scales which only redistributes the kinetic energy between wavelengths. 



III. SKEWNESS OF VELOCITY DIFFERENCE PDF 

n 

The obtained expression of K{r) allows to determine the skewness of Aur ^] 

((Ai,..)^) 6k(r) 

{{Au^m"' {2{i-m)r ' ' 

which is expressed in terms of the longitudinal triple correlation k{r), linked to K{r) by 
K{r) = {d/dr + 4/r) k{r) (also see Appendix A, Eq. P7|) ). Since / and k are, respec- 
tively, even and odd functions of r with /(O) = 1, A;(0) = k'{0) = k"{0) =0, -f^3(0) is given 
by 

i/,(0) = lim/f3(r)=(-^^ (17) 



where the apex denote the derivative with respect to r. To obtain H^lO), observe that, near 
the origin, K behaves as 

K = «3yr7^/"(o)^ + O(r^) (18) 

then, substituting Eq. f|T8l) into K{r) = {d/dr + 4/r) k{r) and accounting for Eq. f|T7j) . 
one obtains 

HsiO) = = -0.42857... (19) 

This value of -^3(0) is a constant of the present analysis, which does not depend on the 
Reynolds number. This is in agreement with the several sources of data existing in the 
literature such as 



0, 13, Q, Q (and Refs. therein) and its value gives the entity of the 



mechanism of energy cascade. 



IV. STATISTICAL ANALYSIS OF VELOCITY DIFFERENCE 

As explained in this section, the Lyapunov analysis of the local deformation and some 
plausible assumptions about the statistics of velocity difference Au(r) = u(X + r) — u(X) 
lead to determine all the statistical moments of Au(r) with only the knowledge of the 
function K{r) and of the value of the critical Reynolds number. 

Starting from the momentum Navier-Stokes equations 

duk duk 1 dTkh 

Ot OXh P OXh 

consider the map x '■ — x, which is the function that__determines the current position x 
of a fluid particle located at the referential position xq 
written in terms of the referential position xq 



at t = to- Equation (120|) can be 



duk _ f duk 1 dTkh \ dxop 



, n + (21) 

ot \ OXop p OXop J OXh 

The Lyapunov analysis of the fluid strain provides the expression of this deformation in 
terms of the maximal Lyapunov exponent 



_ eA(t-to) (22) 

9X0 

where A = A(0) = max(Ai, A2, A3) is the maximal Lyapunov exponent and Aj, {i = 1,2,3) 
are the Lyapunov exponents. Due to the incompressibility, Ai + A2 + A3 = 0, thus, A > 0. 
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If we assume that this deformation is much more rapid than dTkh/dxop and duk/dxopUh, 
the velocity fluctuation can be obtained from Eq. fl2T]) . where dTkh/dxop and duk/dxopUh 
are supposed to be constant with respect to the time 

A V dxop^'" pdxopJt^t^ ^ A\dt ) ^^^^ 

This assumption is justified by the fact that, according to Truesdell jsj], dTkh/dx^p — 
duk/dxopUh is a smooth function of time -at least during the period of a fluctuation- whereas 
the fluid deformation varies very rapidly in proximity of a bifurcation according to Eq. (l22ll . 



The statistical properties of Au(r), are investigated expressing the velocity fluctuation, 
given by Eq. (!23l) . as the Fourier series 

"^lE^w-^-"" (2") 

where U(k) = {Ui{k), U2{i^), U^{k)) are the components of velocity spectrum, which satisfy 
the Fourier transformed Navier-Stokes equations j3] 

^ = -.eupi.)+ 

(25) 

All the components U(k) ~ d\J{K,)/dt/A are random variables distributed according to 
certain distribution functions, which are statistically orthogonal each other 7|. 

Thanks to the local isotropy, u is sum of several dependent random variables which are 
identically distributed [3], therefore u tends to a gaussian variable 16|, and U(k;) satisfies 



the Lindeberg condition, a very general necessary and sufficient condition for satisfying 
the central limit theorem This condition does not apply to the Fourier coefficients of 
Au. In fact, since Au is the difference between two dependent gaussian variables, its PDF 
could be a non gaussian distribution function. In x = 0, the velocity difference Au(r) = 
(Ami, Am2, Ams) is given by 

^"p^X ^ ^^^(e^''" -1)^L + B + P + N (26) 

This fluctuation consists of the contributions appearing into Eq. (125|) : in particular, L 
represents the sum of all linear terms due to the viscosity and B is the sum of all bilinear 



terms arising from inertia and pressure forces. P and are, respectively, the sums of 
definite positive and negative square terms, which derive from inertia and pressure forces. 
The quantity L + B tends to a gaussian random variable being the sum of statistically 
orthogonal terms Q, while P and do not, as they are linear combinations of squares 



171]. Their general expressions are [13] 

P = Po + r]i+ ri 

(27) 

AT = iVo + Ci + CI 

where Pq and Nq are constants, and r/i, r]2i Ci aiid C2 are four different centered random 
gaussian variables. Therefore, the fluctuation Aur of the longitudinal velocity difference can 
be written as 

An, = ^i(r)e + U^) ixiv' " 1) " (C' - 1)) (28) 

where ^, 77 and ( are independent centered random variables which have gaussian distribution 
functions with standard deviation equal to the unity. The parameter x is a positive definite 
function of Reynolds number, whereas ipi and ip2 are functions of space coordinates and the 
Reynolds number. 

At the Kolmogorov scale i, the order of magnitude of the velocity fluctuations is ux'^r/i, 
with r = 1/A and uk = J^/^, whereas 1^2 is negligible because is due to the inertia forces: 
this immediately identifies ipi ~ ux'^r/i. 

On the contrary, at the Taylor scale At, ipi is negligible and the order of magnitude of the 
velocity fluctuations is u'^t/Xt, therefore 1^2 ~ u^t/Xt- 
The ratio ip2/'ipi is a function of R\ 



i^{r,Rx) = ——- ^ — 5-;— = W 1= ^{^) (29) 

^i(r) uk^Xt Y 15715 

where ?/'(r) = 0(1), is a function which has to be determined. 

Hence, the dimensionless longitudinal velocity difference Aur, is written as 

AUr _e + ^(x(^'-l)-(C'-l)) 



(30) 



v/((A^) Vl + 2^2(1 + X') 

The dimensionless statistical moments of Aur are easily calculated considering that ^, t] and 
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( are independent gaussian variables 

^ _ ((A«.)") 



I — n V / 



(31) 



where 



(32) 



1=0 ^ ^ 

In particular, the third moment or skewness, H^, which is responsible for the energy cascade, 
is 

H,^^tit^^ (33) 
(1 + 2'*' (1 + X-')P 

For X 7^ 1) the skewness and all the odd order moments are different from zero, and for n > 3, 
all the absolute moments are rising functions of Rx, thus Aur exhibits an intermittency whose 
entity increases with the Reynolds number. 

All the statistical moments can be calculated once the function x{Rx) ^iid the value of ipQ 
are known. The expression of K{r) obtained in the first part of the work allows to identify 
H^lO) and then fixes the relationship between ipo and x(-Ra) 

_^,(0) = ^*£!(l^ll_ = | (34) 
{1 + 2i,o' {1 + x'd ^ 

where ipo = 4'{0,Rx) = 0{\/Rx) and x = x{R\) > 0. This relationship does not admit 
solutions with x > below a minimum value of {R\)min dependent on ipQ. According to 
the analysis of section HXl (Appendix B), {R\)min is chosen to 10.12, which corresponds to 
'ipo ~ 1.075. (setting x = 0, -Ra = 10.12 in H^{{^)). Varying the value of {R\)min from 8.5 
to 15 would bring values of ipQ between 1.2 and 0.9, respectively. In figure [2], the function 



11 



"I Q I \ I \ I I I \ \ I 

10° 10' 10' 10' 10^ 

FIG. 2: Parameter x plotted as the function of Rx. 

x{R\) is shown for ipQ = 1.075. The hmit x — 0.86592 for i?A C)0 is reached independently 
of the value of ipo- 

The PDF of Aur is expressed through the Frobenious- Perron equation 

F(AO= / / [ piOpiv)p{C) H^<- ^Uri^.vX))d^ dr] dC (35) 

JS, Jri Jc 

where Aur is calculated with Eq. fl30l) . 5 is the Dirac delta and p is a gaussian PDF whose 
average value and standard deviation are equal to and 1, respectively. 

For non-isotropic turbulence or in more complex cases with boundary conditions, the 
velocity spectrum could not satisfy the Lindeberg condition, thus the velocity will be not 
distrubuted following a Gaussian PDF, and Eq. ( l28l) changes its analytical form and can in- 
corporate more intermittant terms [16] which give the deviation with respect to the isotropic 
turbulence. Hence, the absolute statistical moments of Aur will be greater than those cal- 
culated with Eq. fl30l) . indicating that, in a more complex situation than the isotropic 
turbulence, the intermittency of Aur can be significantly stronger. 

V. RESULTS AND DISCUSSION 

In order to obtain informations about the validity of the proposed analysis, several results 
are now presented. 
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As the first result, the evolution in time of the correlation function is calculated with 
the proposed closure of the von Karman-Howarth equation (Eq. f|T5l) ). where the boundary 
conditions are given by Eq. fl49p . The turbulent kinetic energy and the spectrums E{k,) and 
T(k) are calculated with Eq. fISUl) and Eqs. fIFT]) . respectively. The calculation is carried 
out for the initial Reynolds number of Re = u{0)Lr/i' = 2000, where and u{0) are, 
respectively, the characteristic dimension of the problem and the initial velocity standard 
deviation. The initial condition for the correlation function is /(r) = exp (— (r/A^)^), where 
Xx/Lr = l/(2-\/2), whereas ■u(O) = 1. The dimensionless time of the problem is defined as 
t = t u{0)/Lr. 

Equation (1451) was numerically solved adopting the Crank-Nicholson integrator scheme 
with variable time step, where the discretization of the space domain is made by — 1 
intervals of the same amplitude Ar. This corresponds to a discretization of the Fourier 
space made by A^ — 1 subsets in the interval [0, km], where = tt/ (2Ar). For the adopted 
initial Reynolds number, the choice A^ = 1500, gives an adequate discretization, which 
provides Ar < i, for the whole simulation. For what concerns u, it was calculated with Eq. 
( l50l) and the kinetic energy was checked to be equal to the integral over k of the energy 
spectrum. During the simulation, T(k) must identically satisfy Eq. fl52|) (see Appendix A) 
which states that T(k) does not modify the kinetic energy. According to the discretization 
of the Fourier space, the integral of T(/t) is calculated with the trapezes rule from until to 
hm, at each time step, therefore, the simulation will be considered to be accurate as long as 



namely, when T(k) ^ for k > km- As the simulation advances, according to Eq. (|T5|) . 
the energy cascade determines variations of E{k,) and T{k,) for wave-numbers whose values 
rise with the time, then Eq. fl5B]) holds until to a certain time, where these wave-number 
are about equal to km- At higher times, the variations of T(/t) can occur for k > km, out 
of the interval (0, km), thus Eq. (!36l) could be not satisfied. For this reason, the simulation 
is stopped as soon as the following condition is achieved [3] 



At the end of several simulations, we obtain Ar ^ 0.8 i, and, in this situation, the energy 
spectrum is here supposed to be fully developed. 




(36) 




(37) 
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FIG. 3: Correlation functions, / and k versus the separation distance at the times of simulation i 
= 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.63. 

The diagrams of Fig. [3] show the correlation functions /(r) and k{r) vs. the dimensionless 
distance t/Xt, at different times of simulation. The kinetic energy and Taylor scale vary 
according to Eqs. (fT5|) and (!50|) . thus /(r) and k{r) change in such a way that the length 
scales associated to their variations diminish as the time increases, whereas the maximum 
of |/c| decreases. At the final instants of the simulation, one obtains that / — 1 = 0( r^/^) 
for r/A-r = 0(1), whereas the maximum of \k\ is about 0.05. These results are in very 
good agreement with the numerous data of the literature [tI which concern the evolution of 
correlation functions. Figure H] shows the diagrams of E{k,) and T(/t) for the same times, 
where the dashed line in the plot of E{k), represents the —5/3 Kolmogorov law The 
spectrums E{k,) and T(k) vary with time according to Eqs. fll5p and (1511) and depend on the 
initial condition. At the end of simulation, the energy spectrum E{k,) can be compared with 
the dashed line in an opportune interval of wave-numbers. This arises from the developed 
correlation function, which behaves like / — 1 = O (r^/^) for r = 0{\t)- 
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FIG. 4: Plot of E{k) and T{k) at the diverse times of simulation. 
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FIG. 6: (a) Maximum finite size Lyapunov exponent at the times of simulation t = 0, 0.1, 0.2, 0.3, 
0.4, 0.5, 0.6, 0.63; (b) and (c) skewness and Flatness versus r/Ar at t = and t = 0.6, respectively. 

Next, the Kolmogorov function Q{r) and Kolmogorov constant C, are determined with 
the proposed analysis, using the previous results of the simulation. 

Following the Kolmogorov theory, the Kolmogorov function, which is defined as 



Q{r) 



re 



(38) 



is constant with respect to r, and is equal to 4/5 as long as r/Ay = 0(1). As shown in Fig. 
m for t = 0, the maximum of Q{r) is much greater than 4/5 and its variations with r/Xj- can 
not be neglected. This is the consequence of the choice of the initial correlation function. At 
the successive times, the maximum of Q{r) decreases until to the final instants, where, with 
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the exception of r/Xx ~ 0, Q{r) exhibits variations which are less than those calculated at 
the previous times in a wide range of r/Ar, with a maximum which can be compared to 0.8. 
The Kolmogorov constant C is also calculated by definition 

2/3 

E{^) = (39) 

This is here determined, as the value of C which makes the curve represented by Eq. fl39|) to 
be tangent to the energy spectrum E{k) previously calculated. At end simulation, C ~ 1.932, 
namely C and Qmax agree with the corresponding quantities known from the literature. 

For the same simulation. Fig. [6^ shows the maximal finite scale Lyapunov exponent, 
calculated with Eq. (fl^ . where A varies according to /. For t = 0, the variations of A are 
the result of the adopted initial correlation function which is a gaussian, whereas as the time 
increases, the variations of / determine sizable increments of A and of its slope in proximity 
of the origin. Then, for developed spectrum, since / — 1 = 0(r^/^), the maximal finite scale 
Lyapunov exponent behaves like A ~ 7--2/3_ xhus, the diffusivity coefficient associated to 
the relative motion between two fluid particles, defined as D{r) oc Ar^, here satisfies the 
famous Richardson scaling law D{r) ^ r^/^^|. 

In the diagrams of Figs. [6]d and [6t, skewness and flatness of Am^ are shown in terms 
of r for t = and 0.6. The skewness, is first calculated with Eq. ( fT6l) . then has 
been determined using Eq. ( pTl) . At t = 0, l-f^sl starts from 3/7 at the origin with small 
slope, then decreases until to reach small values. also exhibits small derivatives near the 
origin, where ^ 3, thereafter it decreases more rapidly than \H^\. At i =0.6, the diagram 
importantly changes and exhibits different shapes. The Taylor scale and the corresponding 
Reynolds number are both changed, so that the variations of and are associated 
to smaller distances, whereas the flatness at the origin is slightly less than that at t = 0. 
Nevertheless, these variations correspond to higher r/Ar than those for t = 0, and also in 
this case, reaches the value of 3 more rapidly than tends to zero. 

The PDFs of Am^ are calculated with Eqs. (!35l) and (!30l) . and are shown in Fig. [7] in 
terms of the dimensionless abscissa 

Am., 



((An,)2)i/2 

where, these distribution functions are normalized, in order that their standard deviations 
are equal to the unity. The figure represents the distribution functions of s for several r/A^, 
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FIG. 7: PDF of the velocity difference fluctuations at the times t =0 (a), t = 0.5 (b) and t =0.6 
(c). Continuous lines are for r =0, dashed lines are for r/Xx =1, dot-dashed lines are for r/Ar =5, 
dotted lines are for gaussian PDF. 

at t = 0, 0.5 and 0.6, where the dotted curves represent the gaussian distribution functions. 
The calculation of H^{r) is first carried out with Eq. (fT6!) . then the function ip{r,Rx) is 
identified through Eq. fl33l) . and finally the PDF is obtained with Eq. fl35l) . For t = (see 
Fig. and according to the evolutions of and H4, the PDFs calculated at t/At = 
and 1, are quite similar each other, whereas for t/Xt = 5, the PDF is almost a gaussian 
function. Toward the end of the simulation, (see Fig. [7b and c), the two PDFs calculated 
at r/Ar = and 1, exhibit more sizable differences, whereas for r/Ar = 5, the PDF differs 
very much from a gaussian PDF. This is in line with the plots of H^lr) and H^^ir) of Fig. 
IS 

Next, the spatial structure of Am,., given by Eq. fl30l) . is analyzed using the previous 
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15 




FIG. 8: Scaling exponents of longitudinal velocity difference versus the order moment at different 



times. Continuous 
mogorov K41 data_ 
She-Leveque data 



ines with solid symbols are for the present data. Dashed lines are for Kol- 
5]. Dashdotted lines are for Kolmogorov K62 data [l9^. Dotted lines are for 



results of the simulation. According to the various works [l9|, l20|, l2l|, Aur behaves quite 
similarly to a multifractal system, where Aur obeys to a law of the kind Aur{r) ~ where 
the exponent g is a fluctuating function of space. This implies that the statistical moments 
of Aur{r) are expressed through different scaling exponents ({P) whose values depend on 
the moment order P, i.e. 



{iAurfir)) = Ap r«^) 



(40) 



These scaling exponents are here identified through a best fitting procedure, in the intervals 
[ap, ap +Xt), where the endpoints ap are unknown quantities which have to be determined. 
The location of these intervals depends on P and varies with the time. The calculation of 
the endpoints ap and of (p and is carried out through a minimum square method which 
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for each moment order is applied to the following optimization problem 

Jp{apXp,Ap)= / (((A^z,,)0-^pr^(^^)'t^^ = min, P = l,2,... (41) 

J ap 

where ((Am^)) are calculated with Eqs. (PTjl . 

Figure [8] shows the comparison between the scaling exponents here obtained (continuous 
lines with solid symbols) and those of the Kolmogorov theories K41 (dashed lines) and 



K62 [19| (dashdotted lines), and those given by She-Leveque [20| (dotted curves). At t = 
0, the values of C{P) ^'^'^ the result of the chosen initial condition. As the time increases, 
the correlation function changes causing variations in the statistical moments of AMr(r). As 
result, C{P) gradually diminish and exhibit a variable slope which depends on the moment 
order P, until to reach the situation of Fig. [8)d, where the simulation is just ended. The 
dimensionless moments of A-u,.(r) are changed. The plot of C(-P) shows that near the origin, 
C(P) ^ P/3, and that the values of C,{P) seem to be in agreement with the those proposed 
by She-Leveque. More in detail. Table I reports these scaling exponents in terms of the 
moments order, calculated for t = 0.63. These values are the consequence of the spatial 
variations of the skewness, calculated using Eq. (fT6ll . and of the quadratic terms due to 
the inertia and pressure forces into the expression of the velocity difference, which make 
{{Aur)^) a quantity quite similar to a multifractal system. 

Other simulations with different initial correlation functions and Reynolds numbers have 
been carried out, and all of them lead to analogous results, in the sense that, at the end of the 
simulations, the diverse quantities such as Q{r), C and C{P) are quite similar to those just 
calculated. For what concerns the effect of the Reynolds number, its increment determines 
a wider range of the wave-numbers where E{k) is comparable with the Kolmogorov law and 
a smaller dissipation energy rate in accordance to Eq. (l50l) . 

In order to study the evolution of the intermittency vs. the Reynolds number. Table 
II gives the first ten statistical moments of F{dur/dr). These are calculated with Eqs. 
( I3T1) and (l32l) . for Pa = 10.12, 100 and 1000, and are shown in comparison with those of 

P 1 2 3 4 5 6 7 8 9 10 11 12 13 14 

C(P) 0.36 0.71 0.99 1.19 1.41 1.61 1.84 2.04 2.25 2.49 2.72 2.93 3.15 3.38 

TABLE I: Scaling exponents of the longitudinal velocity difference. 
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a gaussian distribution function. It is apparent that a constant nonzero skewness of the 
longitudinal velocity derivative, causes an intermittency which rises with Rx (see Eq. fl30l) ). 
More specifically, Fig. [H] shows the variations of -^4(0) and -^6(0) (continuous lines) in terms 
of R\, calculated with Eqs. flHTl) and fl5^ . with H^lO) = —3/7. These moments are rising 
functions of Rx for 10 < i?A ^ 700, whereas for higher Rx these tend to the saturation 
and such behavior also happens for the other absolute moments. According to Eq. (I3T1) . in 
the interval 10 < i?A ^ 70, H4 and Hq result to be about proportional to R^^"^ and R^'^^, 
respectively, and the intermittency increases with the Reynolds number until to Rx ~ 700, 
where it ceases to rise so quickly. This behavior, represented by the continuous lines, 
depends on the fact that ib ^ ^/Rx, and results to be in very good agreement with the data 
of Pullin and Saffman [221], for 10 < ^ 100. Figure [9] can be compared with the data 
collected by Sreenivasan and Antonia [l5l], which are here reported into Fig. [TOl These latter 
are referred to several measurements and simulations obtained in different situations which 
can be very far from the isotropy and homogeneity conditions. Nevertheless a comparison 
between the present results and those of Ref. 15| is an opportunity to state if the two 



data exhibit elements in common. According to Ref. [15[, the flatness monotonically rises 
with Rx with a rising rate which agrees with Eq. fl32|) for 10 < -Ra ^ 60 (dashed line. 
Fig. [9]), whereas the skewness seems to exhibit minor variations. Thereafter, H4 continues 



Moment Rx ~ 10 Rx = 10^ ^a = 10^ Gaussian 



Order 


P. R. 


P. R. 


P. R. 


Moment 


3 


-.428571 


-.428571 


-.428571 





4 


3.96973 


7.69530 


8.95525 


3 


5 


-7.21043 


-11.7922 


-12.7656 





6 


42.4092 


173.992 


228.486 


15 


7 


-170.850 


-551.972 


-667.237 





8 


1035.22 


7968.33 


11648.2 


105 


9 


-6329.64 


-41477.9 


-56151.4 





10 


45632.5 


617583. 


997938. 


945 



TABLE II: Dimensionless statistical moments of F{dur/dr) at different Taylor scale Reynolds 
numbers. P.R. as for "present results". 
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FIG. 9: Dimensionless moments -^4(0) and -^6(0) plotted vs. R\. Continuous lines are for the 
present results. The dashed line is the tangent to the curve of i^4(0) in R\ ~ 10. 
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FIG. 10: Flatness -^4(0) vs. R\. These data are from Ref. 

to rise with about the same rate, without the saturation observed in Fig. [91 The weaker 
intermittency calculated with the present analysis arise from the isotropy which makes the 
velocity fluctuation a gaussian random variable, while, as seen in sec. IIVI without the 
isotropy condition, the flatness of velocity and of velocity difference can be much greater 
than that of the isotropic case. 



Again, the obtained results are compared with the data of Tabeling et al 13|, |l4|], where, 
in an experiment using low temperature helium gas between two counter-rotating cylinders 
(closed cell), the authors measure the PDF of dur/dr and its moments. Also in this case the 
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FIG. 11: Skewness S = -^3(0), Flatness F = -^4(0) and hyperflatness Hq{0) vs. Rx. These data 
are from Ref . \l\ . 

flow can be quite far from to the isotropy condition. In fact, these experiments pertain wall- 
bounded flows, where the walls could importantly influence the fluid velocity in proximity of 
the probe. The authors found that the higher moments than the third order, first increase 
with Rx until to Rx ~ 700, then exhibit a lightly non-monotonic evolution with respect to 
Rx, and finally cease their variations denoting a transition behavior (See Fig. [TTl) . As far as 
the skewness is concerned, the authors observe small percentage variations. Although the 
isotropy does not describe the non-monotonic evolution near Rx = 700, the results obtained 

n n 

with Eq. (1501) can be considered comparable with those of Refs. [13|, [IJ], resulting also 
in this case, that the proposed analysis gives a weaker intermittency with respect to Refs. 



Il3|, 



1^. 



The normalized PDFs of dur/dr are calculated with Eqs. (!35|) and (!30!) . and are shown 
in Fig. [12] in terms of the variable s, which is defined as 

dur/dr 
{{dur/drff'^ 

Figure [T2k shows the diagrams for Rx= 15, 30 and 60, where the PDFs vary in such a way 
that if3(0) = -3/7. 

As well as in Ref. Q, Figs. 4b and 4c give the PDF for Rx = 255, 416, 514, 1035 and 1553, 
where these last Reynolds numbers are calculated through the Kolmogorov function given 
in Ref. [l^, with H^{0) = —3/7. In particular. Fig. [T2t represents the enlarged region of 
Fig. [T2b . where the tails of PDF are shown for 5 < s < 8. According to Eq. (!30|) . the tails 
of the PDF rise in the interval 10 < -Ra ^ 700, whereas at higher Rx, smaller variations 
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FIG. 12: Log linear plot of the PDF of du^jdr for different R\. (a): dotted, dashdotted and 
continuous lines are for R\ = 15, 30 and 60, respectively, (b) and (c) PDFs for R\ = 255, 416, 
514, 1035 and 1553. (c) represents an enlarged part of the diagram (b) 

occur. Although the non-monotonic trend observed in Ref. IJ], Fig. [T2b shows that the 



values of the PDFs calculated with the proposed analysis, for 5 < s < 8, exhibit the same 
order of magnitude of those obtained by Tabeling et al [l^ which are here shown in Fig. [T51 

Asymmetry and intermittency of the distribution functions are also represented through 
the integrand function of the 4*^* order moment of PDF, which is 0/4(5) = s^F(s). This 
function is shown in terms of s, in Fig. [HI for R\ = 15, 30 and 60. 

VI. CONCLUSIONS 



The proposed analysis is based on the conjecture which states that the turbulence is 
caused by the bifurcations of the velocity field. The main limitation of this analysis is 
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FIG. 13: PDF of dur/dr for = 255, 416, 514, 1035 and 1553. These data are from Ref. [l 
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FIG. 14: Plot of the integrand s'^F{s) for different R\. Dotted, dashdotted and continuous lines 
are for Rx = 15, 30 and 60, respectively. 

that it only studies the developed homogeneous-isotropic turbulence, whereas this does not 
consider the intermediate stages of the turbulence. 
The results of this analysis can be here summarized: 

1. The Lyapunov analysis of the relative motion provides an explanation of the physical 
mechanism of the energy cascade in turbulence and gives a closure of the von Karman- 
Howarth equation. 

The fluid incompressibility determines that the inertia forces transfer the kinetic energy 
between the length scales without changing the total kinetic energy. This implies that 
the skewness of the longitudinal velocity derivative is a constant of the present analysis 
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and that the energy cascade mechanism does not depend on the Reynolds number. 

2. The momentum equations written using the referential coordinates allow to factorize 
the velocity fluctuation and to express it in Lyapunov exponential form of the local 
fluid deformation. The statistics of Aii^ can be inferred looking at the Fourier series of 
the velocity difference. This is a non-Gaussian statistics, where the constant skewness 
of durjdr implies that the other higher absolute moments increase with the Taylor- 
scale Reynolds number. 

3. The momentum equations written using the referential coordinates allow the velocity 
fluctuation to be expressed by means of the Lyapunov analysis of the kinematics of 
local fluid deformation. The Fourier series of the velocity difference provides the 
statistics of A-u^- This is a non-Gaussian statistics, where the constant skewness of 
dur/ dr implies that the other higher absolute moments increase with the Taylor-scale 
Reynolds number. 

4. The closure of the von Karman-Howarth equation, shows that the mechanism of energy 
cascade gives energy spectrums that can be compared with the Kolmogorov law k~^/^ 
in an opportune range of wave-numbers. 

5. For developed energy spectrums, the Kolmogorov function exhibits, in an opportune 
range of r, small variations much less than at the previous times, and its maximum is 
quite close to 4/5, whereas the Kolmogorov constant is about equal to 1.93. As the 
consequence, the maximal flnite scale Lyapunov exponent and the diffusivity coefficient 

vary according to the Richardson law when the separation distance is of the order of 
the Taylor scale. 

6. The analysis also determines the scaling exponents of the moments of the longitudinal 

velocity difference through a best fitting procedure. For developed energy spectrum, 
these exponents show variations with the moment order which seem to be consistent 
with those known from the literature. 
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VIII. APPENDIX A 

The von Karman-Howarth equation gives the evolution in time of the longitudinal corre- 
lation function for isotropic turbulence. The correlation function of the velocity components 
is the symmetrical second order tensor Rij{r) = iuiu'^^^ where Ui and u'j are the velocity 
components at x and x + r, respectively, being r the separation vector. The equations for 



Rij are obtained by the Navier-Stokes equations written in the two points x and x + r 
For isotropic turbulence Rij can be expressed as 



Rijir) = 



if-9f-^ + 9S^, (42) 



/ and g are, respectively, longitudinal and lateral correlation functions, which are 
^^^^ ^ {ur{x)ur{x + r)) ^^^^ ^ (n„ (r)M„(x + r) ) 

where Ur and m„ are, respectively, the velocity components parallel and normal to r, whereas 
r = |r| and = {uf) ={uD= 1/3 {utUi). Due to the continuity equation, / and g are linked 
each other by the relationship 

9 = f^l%r (44) 



The von Karman-Howarth equation reads as follows 



m 



df K fd^f 4df\ d^f , 

lt = - + ^^{^ + Zlt]- 10^7^(0)/ (45) 



dt \ dr"^ r dr J dr 

where K is an even function of r, which is defined by the following equation 



m 



and which can also be expressed as 



Kir) = «M |: + ^ ) Hr) (47) 
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where k is the longitudinal triple correlation function 

(M^(x)ur(x + r)) 



k{r) 



m 



(48) 



The boundary conditions of Eq. fH5l) are 

/(0) = 1, lim/(r) = 



(49) 



The viscosity is responsible for the decay of the turbulent kinetic energy, the rate of which 

is BEi 



(50) 



This energy is distributed at different wave-lengths according to the energy spectrum E{k) 
which is calculated as the Fourier Transform of fv^, whereas the "transfer function" Tin) 

n 

is the Fourier Transform of i^T [7|, i.e. 



E{k) 








1 




Trio 




T{k) _ 




K{r) 



2 2 / sm Kr 
K r I cos Kr I dr 



(51) 



where k = |k;| and T{k,) identically satisfies to the integral condition 

T{K)dK = (52) 

which states that K does not modify the total kinetic energy. The rate of energy dissipation 
e is calculated for isotropic turbulence as follows 



3du^ P 2r./ X , 

e = — =2v I K EiK)dK 

2 dt Jo ^ ^ 

The microscales of Taylor At, and of Kolmogorov i, are defined as 



(53) 



A 



u 



3\ 1/4 



(54) 
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IX. APPENDIX B: CRITICAL REYNOLDS NUMBER 



The purpose of this appendix is to provide an estimation of the critical Reynolds number 
assuming that the turbulence is, in any case, fully developed, homogeneous and isotropic. 
Thus, the obtained results are subjected to these assumptions. 

To this end, consider now the equation of motion of a fluid particle c/x/c/t = u(x, t) and 
its fixed points which satisfy rfx/rft = 0. We assume that the bifurcations cascade of this 
equation are expressed in terms of the characteristic scales by the asymptotic approximation 



(55) 



where a ~ 2 [8|, l23|, and represent the average distance between two branches of fixed 
points which born in the same bifurcation. Equation (1551) is supposed to describe the route 
toward the chaos and is assumed to be valid until the onset of the turbulence. In this 
situation the minimum for /„ can not be less than the dissipation length or Kolmogorov 
scale ^ = (i^^/fV^ [H, where h gives a good estimation of the correlation length of the 
phenomenon [sl, which, in this case is the Taylor scale Xt- Thus, i < In < \t, and 

" (56) 



where is the number of bifurcations at the beginning of the turbulence. Equation ( l56l) gives 
the connection between the critical Reynolds number and number of bifurcations. In fact, the 
characteristic Reynolds numbers associated to the scales i and Xt are Rk = £uk/i^ = 1 and 
R\ = Xtu/u, respectively, where uk = (z^e)^''^ is characteristic velocity at the Kolmogorov 
scale, and u = \J (uiUi) /3 is the velocity standard deviation [7]. For isotropic turbulence, 
these scales are linked each other by [3] 

Ay/£= (57) 

In view of Eq. (1561) . this ratio can be also expressed through N, i.e. 



-1 = 151/^/Rl (58) 



Assuming that a is equal to the Feigenbaum constant (2.502...), the value Rx ~ 1.6 obtained 
for = 2 is not compatible with which is the correlation scale, while the result Rx ~ 

29 



10.12, calculated for = 3, is an acceptable minimum value for Rx. The order of magnitude 
of these values can be considered in agreement with the various scenario describing the roads 
to the turbulence 
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26l |. and with the diverse experiments 
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29( 1 which state 



that the turbulence begins for > 3. Of course, this minimum value for Rx is the result of 



the assumptions a ~ 2.502, li ~ At, In — ^ and of the asymptotic approximation (ITOI) . 
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